####Table A9####
main_m1dd_t1_sb <- glm(as.formula(paste("osv_event", paste(c(m1d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","osv_event_peaceyears_l1","I(osv_event_peaceyears_l1^2)","I(osv_event_peaceyears_l1^3)","osv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), data=subset(main_ddyad, year <= 2012 & osv_sample == 1 & int_grp_rel_dominant_g1 == 1), family=binomial(link="logit"),control = list(maxit = 50))
main_m1dd_t1_sb_cse <- cluster.se(main_m1dd_t1_sb, as.factor(subset(main_ddyad, year <= 2012 & osv_sample == 1 & int_grp_rel_dominant_g1 == 1)$cowcode))
main_m2dd_t1_sb <- glm(as.formula(paste("osv_event", paste(c(m2d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","osv_event_peaceyears_l1","I(osv_event_peaceyears_l1^2)","I(osv_event_peaceyears_l1^3)","osv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), data=subset(main_ddyad, year <= 2012 & osv_sample == 1 & int_grp_rel_dominant_g1 == 1), family=binomial(link="logit"),control = list(maxit = 50))
main_m2dd_t1_sb_cse <- cluster.se(main_m2dd_t1_sb, as.factor(subset(main_ddyad, year <= 2012 & osv_sample == 1 & int_grp_rel_dominant_g1 == 1)$cowcode))
main_m3dd_t1_sb <- glm(as.formula(paste("osv_event", paste(c(m3d_t1, dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","osv_event_peaceyears_l1","I(osv_event_peaceyears_l1^2)","I(osv_event_peaceyears_l1^3)","osv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), data=subset(main_ddyad, year <= 2012 & osv_sample == 1 & int_grp_rel_dominant_g1 == 1), family=binomial(link="logit"),control = list(maxit = 50))
main_m3dd_t1_sb_cse <- cluster.se(main_m3dd_t1_sb, as.factor(subset(main_ddyad, year <= 2012 & osv_sample == 1 & int_grp_rel_dominant_g1 == 1)$cowcode))
stargazer(main_m1dd_t1_sb, main_m2dd_t1_sb, main_m3dd_t1_sb, se=c(data.frame(main_m1dd_t1_sb_cse[, 2]),data.frame(main_m2dd_t1_sb_cse[, 2]), data.frame(main_m3dd_t1_sb_cse[, 2])), single.row = T, dep.var.labels.include = F, omit=c("cowcode|region|factor"), type="text", model.numbers =F, column.labels = c("model A162","model A163","model A164"), order=vars.order_d, title="Table A9. Territorial autonomy and one-sided violence incidence by second-order majorities against second-order minorities.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of one-sided violence perpetrated by the given second-order majority against the given second-order minority in a given unit.", "Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Territorial autonomy x included (maj.)/excluded (min.)","Territorial autonomy x included (min.)/excluded (maj.)","Included/excluded","Included (maj.)/excluded (min.)","Included (min.)/excluded (maj.)","Excluded/excluded","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","One-sided violence peace years", "One-sided violence peace years 2", "One-sided violence peace years 3","One-sided violence spatial lag"), out="../tables/tablea9.txt")
stargazer(main_m1dd_t1_sb, main_m2dd_t1_sb, main_m3dd_t1_sb, se=c(data.frame(main_m1dd_t1_sb_cse[, 2]),data.frame(main_m2dd_t1_sb_cse[, 2]), data.frame(main_m3dd_t1_sb_cse[, 2])), single.row = T, dep.var.labels.include = F, omit=c("cowcode|region|factor"), type="html", model.numbers =F, column.labels = c("model A162","model A163","model A164"), order=vars.order_d, title="Table A9. Territorial autonomy and one-sided violence incidence by second-order majorities against second-order minorities.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of one-sided violence perpetrated by the given second-order majority against the given second-order minority in a given unit.", "Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Territorial autonomy x included (maj.)/excluded (min.)","Territorial autonomy x included (min.)/excluded (maj.)","Included/excluded","Included (maj.)/excluded (min.)","Included (min.)/excluded (maj.)","Excluded/excluded","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","One-sided violence peace years", "One-sided violence peace years 2", "One-sided violence peace years 3","One-sided violence spatial lag"), out="../tables/tablea9.html")

####Figure 6b####
#note: figure 6 is the combination of figure 6, panels a, b, and c
me_main_m1dd_t1_sb <- margins_summary(main_m1dd_t1_sb, variables = "sa_territory_t", vcov = cluster.vcov(main_m1dd_t1_sb, as.integer(subset(main_ddyad, year <= 2012& int_grp_rel_dominant_g1 == 1 & osv_sample == 1)$cowcode)))
me_main_m1dd_t1_sb$term <- "all"
me_main_m1dd_t1_sb$model <- 1
me_main_m2dd_t1_sb <- margins_summary(main_m2dd_t1_sb, variables = "sa_territory_t", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(main_m2dd_t1_sb, as.integer(subset(main_ddyad, year <= 2012& int_grp_rel_dominant_g1 == 1 & osv_sample == 1)$cowcode)))
me_main_m2dd_t1_sb$term <- ifelse(me_main_m2dd_t1_sb$included_excluded == 1, "included/excluded", "other")
me_main_m2dd_t1_sb$model <- 1
me_main_m3dd_t1_sb1 <- margins_summary(main_m3dd_t1_sb, variables = "sa_territory_t", at = list(includedd_excludednd = 1, includednd_excludedd = 0), vcov = cluster.vcov(main_m3dd_t1_sb, as.integer(subset(main_ddyad, year <= 2012& int_grp_rel_dominant_g1 == 1 & osv_sample == 1)$cowcode)))
me_main_m3dd_t1_sb1$term <- "included (maj.)/excluded (min.)"
me_main_m3dd_t1_sb1$model <- 1
me_main_m3dd_t1_sb2 <- margins_summary(main_m3dd_t1_sb, variables = "sa_territory_t", at = list(includednd_excludedd = 1, includedd_excludednd = 0), vcov = cluster.vcov(main_m3dd_t1_sb, as.integer(subset(main_ddyad, year <= 2012& int_grp_rel_dominant_g1 == 1 & osv_sample == 1)$cowcode)))
me_main_m3dd_t1_sb2$term <- "included (min.)/excluded (maj.)"
me_main_m3dd_t1_sb2$model <- 1
me_main_m13dd <- rbind.fill(me_main_m1dd_t1_sb, me_main_m2dd_t1_sb, me_main_m3dd_t1_sb1,me_main_m3dd_t1_sb2)
me_main_m13dd <- subset(me_main_m13dd, term != "other")
me_main_m13dd$term <- as.factor(me_main_m13dd$term)
me_main_m13dd$term = factor(me_main_m13dd$term,levels(me_main_m13dd$term)[c(1,3,4,2)])
me_main_m13dd$estimate <- me_main_m13dd$AME
me_main_m13dd$conf.low <- me_main_m13dd$lower
me_main_m13dd$conf.high <- me_main_m13dd$upper
figure6b <- dwplot(me_main_m13dd,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(shape = 16, colour = "black"), whisker_args = list(linetype = 1, colour = "black")) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("b) AME of territorial autonomy on one-sided violence maj./min. dyad") +
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of one-sided violence") + ylab("dyad type") +
  scale_x_continuous(labels=scales::percent_format())
ggsave(figure6b, file='../figures/figure6b.pdf', width = 14, height = 3.5, units="cm",dpi=1000)